Multi-reference symmetry-projected variational approaches for ground and excited 

states of the one-dimensional Hubbard model 
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We present a multi-reference configuration mixing scheme for describing ground and excited states, 
with well defined spin and space group symmetry quantum numbers, of the one-dimensional Hubbard 
model with nearest-neighbor hopping and periodic boundary conditions. Within this scheme, each 
state is expanded in terms of non- orthogonal and variationally determined symmetry-projected 
configurations. The results for lattices up to 30 and 50 sites compare well with the exact Lieb- 
Wu solutions as well as with results from other state-of-the-art approximations. In addition to 
spin-spin correlation functions in real space and magnetic structure factors, we present results for 
spectral functions and density of states computed with an ansatz whose quality can be well-controlled 
by the number of symmetry-projected configurations used to approximate the systems with JV e 
and N e ± 1 electrons. The intrinsic symmetry-broken determinants resulting from the variational 
calculations have rich structures in terms of defects that can be regarded as basic units of quantum 
fluctuations. Given the quality of the results here reported, as well as the parallelization properties 
of the considered scheme, we believe that symmetry-projection techniques, which have found ample 
applications in nuclear structure physics, deserve further attention in the study of low-dimensional 
correlated many-electron systems. 

PACS numbers: 71.27.+a, 74.20.Pq, 71.10.Fd 



I. INTRODUCTION. 

Studies of correlations arising from electron-electron 
interactions remain a central theme in condensed mat- 
ter physisi to better understand challenging phenomena 
such as high-Tc superconductivity 2 or colossal magnetic 
resistance^ There is a need for better theoretical models 
that can account for relevant correlations in ground and 
excites states of fermionic systems with as much sim- 
plicity as possible. Within this context, the repulsive 
Hubbard Hamiltonian^ has received a lot of attention 
since it is considered the generic model of strongly cor- 
related electron systems. 1 Hubbard- like models have also 
received renewed attention in the study of cold fermionic 
atoms in optical lattices^ and the electronic properties of 
graphene. 6 

Unlike the one-dimensional (ID) Hubbard model, 
which is exactly solvable 7 - using the Bethe ansatz, 8 an 
exact solution of the two-dimensional (2D) problem is 
not known. Therefore, it is highly desirable to develop 
approximations that, on the one hand, can capture the 
main features of the exact ID Bethe solution and, on the 
other hand, can be extended to higher dimensions. For 
small lattices, one can resort to exact diagonalization us- 
ing the Lanczos methodic For larger systems, several 
other methods have been extensively used to study the 
ID and 2D Hubbard models as well as their strong cou- 
pling versions. 10 Among such approximations, we have 
the quantum Monte Carlo ) 11 : 12 the variational Monte 
Carlo,— the density matrix renormalization groupi^— as 
well as approximations based on matrix product and ten- 
sor network states. 18 Both the dynamical mean field the- 
ory and its cluster extensions^— have made important 



contributions to our present knowledge of the Hubbard 
model. Other embedding approaches are also available. 27 
Finally, we refer the reader to the recent state-of-the-art 
applications of the coupled cluster method to frustrated 
Hubbard-like models i 28 ' 29 

Although routinely used in nuclear structure physics, 
especially within the Generator Coordinate Method ) 30 : 31 
symmetry restoration via projection techniques 3 ^ has re- 
ceived little attention in condensed matter physics. Nev- 
ertheless, these techniques offer an alternative for ob- 
taining accurate correlated wave functions that respect 
the symmetries of the considered many-fermion prob- 
lem. The key idea is^ 2 - on the one hand, to consider 
a mean-field trial state \T>) which deliberately breaks 
several symmetries of the original Hamiltonian. On the 
other hand, the Goldstone manifold R\T>), where R rep- 
resents a symmetry operation, is degenerate and the 
superposition of such Goldstone states^ can be used 
to recover the desired symmetry by means of a self- 
consistent variation-after-projection procedure! 32 ' 33 Such 
a single-reference (SR) scheme provides the optimal Ritz- 
variational— representation of a given state by means 
of only one symmetry-projected mean-field configura- 
tion. This kind of SR variation-after-projection scheme, 
has already been applied to the ID and 2D Hubbard 
model o 35 ' 36 as well as in quantum chemistry within the 
framework of the Projected Quasiparticle Theory^ - — 

One of the main advantages of the symmetry-projected 



approximations 



32.35-39 



is that they offer compact wave 



functions as well as a systematic way to improve their 
quality by adopting a multi-reference (MR) approach. In 
this case, a set of symmetry-broken mean-field states \T> 1 ) 
is used to build Goldstone manifolds R\V l ) whose super- 



position can be used to recover the desired symmetries 
of the Hamiltonian i 40 i 41 The key idea is then to expand 
a given state in terms of several symmetry-projected and 
variationally determined mean-field configurations. The 
resulting wave functions encode more correlations than 
the ones obtained within SR methods while still keeping 
well defined symmetry quantum numbers^ 2 - 

There are differents flavors of MR approximations 
available in the literature42r— In the present study, we 
adopt a MR scheme well known in nuclear structure 
physics, - 4C) which to the best of our knowledge has not 
been applied to lattice models. The key ingredient in 
such MR scheme is the inclusion of relevant correlations 
in both ground and excited states on an equal footing. 
As a benchmark test, we concentrate on the ID Hubbard 
model for which exact solutions are known££ In particu- 
lar, wc consider the case of half-filled lattices. Neverthe- 
less, the present MR approximation can be extended to 
the 2D case as well as to doped systems with arbitrary 
on-site interaction strengths. 

For a given single-electron space we resort to gener- 
alized HF-transformations 50 (GHF) mixing all quantum 
numbers of the single-electron basis states. The corre- 
sponding Slater determinants deliberately break spin and 
spatial symmetries of the ID Hubbard Hamiltonian^ 
We restore these broken symmetries with the help of pro- 
jection operators^ The resulting MR ground state wave 
functions are obtained applying the variational principle 
to the projected energy. 

The structure of our MR ground state wave functions is 
formally similar to the one adopted within the Resonat- 
ing HF— ~— (ResHF) method, i.e., they are expanded in 
terms of a given number of non-orthogonal symmetry- 
projected configurations. Nevertheless, while in the lat- 
ter all the underlying HF-transformations and mixing co- 
efficients are optimized simultaneously) 42 ' 46 in our case 
the orbital optimization is performed sequentially, only 
for the last added HF-transformation (all our mixing co- 
efficients are still optimized at the same time) rendering 
our calculations easier to handle. This is particularly rel- 
evant for alleviating our numerical effort if one keeps in 
mind that, for both ground and excited states, we use 
the most general GHF-transformations and therefore a 
full 3D spin projection is required. 

Our MR scheme is also used to compute spin-spin cor- 
relation functions (SSCFs) in real space, magnetic struc- 
ture factors (MSFs) as well as dynamical properties of 
the ID Hubbard model like spectral functions (SFs) and 
density of states (DOS) Ja25i2SiSl On the other hand, one 
may wonder whether there is any relevant information in 
the intrinsic symmetry-broken GHF-determinants asso- 
ciated with our MR wave functions. As will be shown 
below, the structure of such intrinsic determinants can 
be interpreted in terms of basic units of quantum fluctu- 
ations for the lattices considered. 52 

In addition to ground state properties, our MR 
framework treats excited states with well defined 
quantum numbers as expansions in terms of non- 



orthogonal symmetry-projected configurations using 
chains of variation-after-projection (VAP) calculations. 
As a byproduct, we also obtain a (truncated) basis con- 
sisting of a few Gram-Schmidt orthonormalized states^ 
which may be used to perform a final diagonalization of 
the Hamiltonian in order to account for further correla- 
tions in both ground and excited states. 

The layout of the theory part of this paper is as follows. 
First, we introduce the methodology of our MR VAP 
scheme in SecJlTJ Symmetry restoration based on a sin- 
gle Slater determinant (i.e., SR symmetry restoration) is 
described in Sec lII Al This section will serve to set our no- 
tation as well as to introduce some key elements of our 3D 
spin and full space group projection techniques. Subse- 
quently, symmetry restoration based on several Slater de- 
terminants (i.e., MR symmetry restoration) is discussed 
in Sec lII Bl In particular, the MR description of ground 
and excited states is presented in Secs lIIBTI and III F3 21 
respectively. In Sec lII Cl we will briefly discuss the com- 
putation of the SFs and DOS within our theoretical 
framework. 

The results presented in this paper test the perfor- 
mance of our approximation in a selected set of illus- 
trative examples. In most cases, calculations have been 
carried out for on-site repulsions U = 2t, At and 8i taken 
as representatives of weak, intermediate-to-strong, and 
strong correlation regimes. In Sec lIIII we first consider 
the ground states of half-filled lattices with up to 50 
sites. We compare our ground state and correlation en- 
ergies with the exact ones as well as with those obtained 
using other theoretical methods. We then discuss the 
dependence of the predicted correlation energies on the 
number of non-orthogonal symmetry-projected configu- 
rations used to expand our ground state wave functions, 
the computational performance of our scheme as well as 
the structure of the intrinsic GHF-determinants result- 
ing from our MR VAP procedure. Next, we consider the 
results of our calculations for SSCFs in real space and 
MSFs for half-filled lattices with up to 30 sites. Subse- 
quently, we compare the DOS provided by our theoretical 
framework with the exact one, obtained with an in-house 
diagonalization code, in a lattice with 10 sites. Results 
for hole SFs are also discussed for a 30-site lattice. We 
end Sec lIIII by presenting results for excitation spectra 
in various lattices and discussing the structure of the in- 
trinsic GHF-determinants resulting from our MR VAP 
procedure for excited states. Finally, SeclIVI is devoted 
to concluding remarks and work perspectives. 



II. THEORETICAL FRAMEWORK 

In what follows, we describe the theoretical framework 
used in the present study. First, SR symmetry restora- 
tion is presented in Sec lII Al Subsequently, in Sec lII Bl 
we consider our MR scheme to describe both ground 
( Sec lII Bip and excited (Sec lII B2|) states of the ID Hub- 
bard model. The computation of SFs and DOS is briefly 



discussed in Sec lII CI 

A. Single-reference (SR) symmetry restoration 

We consider the ID Hubbard Hamiltonian 4 
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where T fc , (a) is the matrix representation of an irre- 
ducible representation, which can be found by standard 
methods,—^ 2 - and R(g) represents the corresponding sym- 
metry operations (i.e., translation by one lattice site and 
the reflection x — > —x) parametrized in terms of the la- 
bel g. The linear momentum k — (2it / N S i tes ) £, is given 
in terms of the quantum number £ that takes the values 



N s 



1,. 



N s 



(5) 



where the first term represents the nearest-neighbor hop- 
ping (t > 0) and the second is the repulsive on-site inter- 
action (U > 0). The fermionio^i operators c] CT and Cj a 
create and destroy an electron with spin-projection a = 
±1/2 (also denoted as a —"[, I) along an arbitrary chosen 
quantization axis on a lattice site j = 1, . . . , N S it es - The 
operators n Jcr = c'-Cj a are the local number operators. 
We assume periodic boundary conditions, i.e., the sites j 
and j + N s ites are identical. Furthermore, we assume a 
lattice spacing A = 1. 

In the standard HF-approximation ; 32 ' 34 the ground 
state of an 7V e -electron system is represented by a Slater 
determinant \V) — ]~Ii=i bh\0) m which the energetically 
lowest N e single-fermion states (holes h, h , . . . ) are oc- 
cupied while the remaining 2N S i tes — N e states (particles 
p,p , . . . ) are empty. For a set of single-fermion opera- 
tors & , the HF-quasiparticle operators b' are given by 
the following canonical transformatio n 32 ' 34 



allowed inside the Brillouin zone (BZ)i^ Equivalently, it 
can take all integer values between and N S i tes ~ 1 ■ For 
k = 0, 7r an additional label b = ±1 should be introduced 
to account for the parity of the corresponding irreducible 
representation under the reflection x — > —x£££ In what 
follows, we do not explictly write this label b but the 
reader should keep in mind that it is taken into account 
whenever needed. 

We introduce 3 - 6 - the shorthand notation = (S, k) for 
the set of symmetry (i.e., spin and linear momentum) 
quantum numbers as well as K = (E,m). The total 
projection operator reads 



pG _ pS pk 
1 KK' — r T.Y,' r mm 



(6) 



We then superpose the Goldstone manifold 
R(Q,)R(g)\T>) to recover the spin and spatial 
symmetries 3 - 2 " via the following SR ansatz 



b] 






(2) 



where T> is a general 2N S it es X 2iV S i tes unitary 53 matrix, 
i.e., VD^ = T>*T> — 1. In all the calculations to be 
discussed below, we have used generalized HF (GHF) 
transformations. 50 As it is well known, the most general 
GHF-determinant \T>) deliberately breaks several sym- 
metries of the original Hamiltonia m 32 ' 35-38 ' 40 Typical ex- 
amples are the rotational (in spin space) and spatial 
symmetries. To restore the spin quantum numbers in 
a symmetry-broken GHF-determinant, we explicitly use 
the full 3D projection operator— ~— 
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2S+1 
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dttDli,(£i)R(£i) 
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where R(ft) 



-iaS Z( ,- 



-i/3S 



v e 



-ilS z 



is the rotation opera- 



tor in spin space, the label ft — (a, /3, 7) stands for the 
set of Eulcr angles and 2?™/ (&) are Wigner functions^ 
To recover the spatial symmetries, we introduce the pro- 
jection operator 



\V-Q-K)=Y^f%P% K -\V) 



(7) 



K' 



where / e are variational parameters. Note, that the 
state Eq.([7]) is already multi-determinanta l 36 ' 37 via the 
projection operator 2^!x , ■ For a given symmetry 0, 
the energy (independent of K) associated with the state 
Eq.0 



E { - 
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jetTv-o/e 
is given in terms of the Hamiltonian and norm 
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matrices. It has to be minimized with respect to the 
coefficients / e and the underlying GHF-transformation 
V. The variation with respect to the former yields the 
following resonon-like 5 - 6 - eigenvalue equatio n 35 ' 36 
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TABLE I: Ground state energy of lattices with N sites — 12 
and 20, as predicted with the GHF-FED scheme based on 
n\ — 10 GHF-transformations, are compared with exact re- 
sults for on-site repulsions of U — 2t, At and 8t. Energies ob- 
tained with the RHF and UHF approximations are included 
as a reference. The ratio of correlation energies k obtained 
with the UHF and GHF-FED aproximations, is computed ac- 
cording to Eq. (|27[) . For more details, see the main text. 









Nsites — 12 


K{%) 


Nsites — 20 


K{%) 


u= 


=2t 


RHF 


-8.9282 




-15.2551 








UHF 


-9.3379 


36.79 


-15.6411 


24.05 






GHF-FED 


-10.0401 


99.85 


-16.8565 


99.79 






EXACT 


-10.0418 




-16.8599 




u= 


=4t 


RHF 


-2.9282 




-5.2550 








UHF 


-5.6290 


67.65 


-9.3821 


66.08 






GHF-FED 


-6.9201 


99.99 


-11.4954 


99.92 






EXACT 


-6.9204 




-11.5005 




u= 


=8t 


RHF 


9.0718 




14.7450 








UHF 


-2.9532 


92.26 


-4.9219 


92.23 






GHF-FED 


-3.9625 


99.99 


-6.5612 


99.96 






EXACT 


-3.9626 




-6.5699 





with the constraint / t_/V e y e — \ ensuring the orthog- 
onality of the solutions. On the other hand, the unre- 
stricted minimization of the energy [Eg. ([8])] with respect 
to T> is carried out via the Thouless theoremi 35 ' 36 ' 53 

For a given symmetry O, we only retain the en- 
ergetically lowest solution of our VAP equations. 36 
Both the GHF-transformation T> and the mixing coef- 
ficients / e are complex, therefore one needs to mini- 
mize n var = 2(2N sites — N e ) x N e + AS real variables. 
We use a limited-memory quasi-Newton method for such 
minimization . 35 ' 36 ' 57 In practice, the integration over the 
set of Euler angles in Eq.Q is discretized. For example, 
for a lattice with N S it es — 30 we have used 13, 26, and 13 
grid points for the integrations over a, /3, and 7, respec- 
tively. In this case, a total of 263,640 grid points are used 
in the discretization of the projection operator of Eq.©. 
We have afforded such a task by developing a parallel im- 
plementation for all the VAP schemes discussed in this 
paper. 



known in nuclear structure physics as the FED VAMP 40 
(Few Determinant Variation After Mean-field Projec- 
tion) strategy, for the considered ground states are de- 
scribed in the next subsection. We use the acronym 
GHF-FED to refer to it in the present work. On the other 
hand, our MR approach for excited states, known as EX- 
CITED FED VAMP, 40 will be presented in Sec llTTHl 
We will use the acronym GHF-EXC-FED to refer to it 
in what follows. 



1. MR symmetry restoration for ground states (GHF-FED) 

Our goal in this section is to obtain, through a chain 
of VAP calculations, non-orthogonal symmetry-projected 
GHF-configurations used to build a MR expansion of a 
given ground stated with well defined symmetry quan- 
tum numbers O. 

Suppose we have generated a ground state solution 
|^i|) = \V\;Q;K) [Eq.0]. Note that at this point, 
we have added the superscript 1 to explicitly indicate 
that only one GHF-transformation has been used within 
the SR approximation discussed in Sec lII Al On the other 
hand, the subscript 1 has been added to indicate that the 
ground state is considered. As we will see in Sec lII B 21 
this subscript will allow us to distinguish between ground 
(i.e., i = l) and excited (i.e., i = 2, 3, . . . , m) states. On 
the other hand, both indices are also added to the (in- 
trinsic) GHF-transformation to explicitly indicate that it 
is variationally optimized for the state \4>\% ). We then 
keep the transformation T>\ fixed and consider the ansatz 



K' i=1 



(11) 



which approximates the ground state (subscript 1) by 
means of two (superscript 2) non-orthogonal symmetry- 
projected GHF-determinants. It is obtained applying the 
variational principle to the energy functional with re- 
spect to the last added transformation T>\ and all the 
new mixing coefficients f{ B . A similar procedure can 
be followed to approximate the ground state by a larger 
number of non-orthogonal symmetry-projected configu- 
rations. Let us assume that ni — 1 configurations have al- 
ready been computed. Then, one introduces a new GHF- 
transformation 2?" 1 , a new set of mixing coefficients j\ , 
and makes the MR GHF-FED ansatz 



B. Multi-reference (MR) symmetry restoration 

For each symmetry O, the SR procedure described in 
Sec lII Al provides us with the optimal variational repre- 
sentation of the corresponding ground state via a sin- 
gle symmetry-projected GHF-determinant. However, as 
the lattice size increases one may adopt a MR perspec- 
tive to keep and/or improve the quality of the wave 
functions] 40 ' 41 The key features of our MR approach, 



\4>Tk)=Y,Y,^k' P kk'\ D \) 

k' i=1 



E 



me 



fniGfajraiO fniO 



(12) 



which superposes the Goldstonc manifolds 
R(fl)R(g)\T>\). The corresponding energy 



(13) 



TABLE II: Ground state energy of lattices with N S it es = 30 
and 50 predicted with the GHF-FED scheme based on rii — 
25 GHF-transformations, are compared with exact results for 
on-site repulsions of U — 2t, At and 8t. Results obtained with 
the UHF-ResHF approximation,— based on ni = 30 UHF- 
transformations, as well as the RHF and UHF energies are 
also included in the table. The ratio of correlation energies 
k obtained with the UHF, UHF-ResHF and the GHF-FED 
aproximations, is computed according to Eg. (1271) . 



N a . 



= 30 «;(%) N s 



= 50 



<%) 



U=2t 


RHF 


-23.2671 




-38.7039 






UHF 


-23.4792 


10.02 


-39.1294 


12.02 




UHF-ResHF 


-25.3436 


98.11 


-41.9535 


91.78 




UHF-FED 


-25.3508 


98.45 


-41.9963 


92.99 




GHF-FED 


-25.3730 


99.50 


-42.1219 


96.46 




EXACT 


-25.3835 




-42.2443 




U=4t 


RHF 


-8.2671 




-13.7039 






UHF 


-14.0732 


64.75 


-23.4553 


65.02 




UHF-ResHF 


-17.0542 


98.00 


-27.9633 


95.09 




UHF-FED 


-16.9420 


96.75 


-27.3518 


91.01 




GHF-FED 


-17.1789 


99.39 


-27.9788 


95.19 




EXACT 


-17.2335 




-28.6993 




U=8t 


RHF 


21.7329 




36.2961 






UHF 


-7.8329 


93.65 


-12.3048 


92.26 




UHF-ResHF 


-9.5378 


99.04 


-15.6422 


98.59 




UHF-FED 


-9.3524 


98.46 


-14.8461 


97.08 




GHF-FED 


-9.7612 


99.75 


-15.6753 


98.65 




EXACT 


-9.8387 




-16.3842 





is given in terms of the Hamiltonian and norm 



iK,jK' 
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iKJK 
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kernels, which require the knowledge of the symmetry- 
projected matrix elements between all the GHF- 
determinants used in the expansion Eg. (1121) . The 
wave function Eq. (|12p is determined varying the energy 
Eq. ([T3|) with respect to all the new mixing coefficients 
fl and the last added transformation I?™ 1 . In the for- 
mer case, we obtain an eigenvalue equation similar to 
Eq.flHD, with the constraint /mQtjV"i©/nie _ ^ while 
the unrestricted minimization with respect to Z?™ 1 is car- 
ried out via the Thouless theorem. Let us stress that the 
GHF-FED MR approximation Eq fT2l) of a given ground 
state enlarges the flexibity in our wave functions to a total 
number of n var = 2ni(2N S i tes — N e )x N e +4:niS+2(ni — 1) 
variational parameters. 



2. MR symmetry restoration for excited states 
(GHF-EXC-FED) 

In this section, we construct non-orthogonal 
symmetry-projected GHF-configurations to expand 
a given excited state. The orthogonalization be- 
tween ground and excited states is achieved via the 
Gram-Schmidt procedure. 36 As a byproduct, our MR 
GHF-EXC-FED method also yields a (truncated) basis 
consisting of a few orthonormal states which may be used 
to diagonalize the Hamiltonian and account for further 
correlations in both ground and excited states^ 6 -^ 

Let us assume that we have already obtained a GHF- 
FED ground state I0'/ 1 } = |<ffl^ e ) [Eq HH)|] along the 
lines discussed in the previous Sec lIIB II We then look 
for the first excited state (subscript 2) with the same 
symmetry 0, approximated by a given n^ number of non- 
orthogonal symmetry-projected configurations. We start 
with the ansatz 



\<d)=o?\W)+P 1 \<&) 



(15) 



where \<f>2%) has a form similar to Eq.([7|) but written in 
terms of the coefficients j\® and the GHF-determinant 
\T>2). Both a 1 and [3 1 can be obtained by requiring or- 
thonormalization with respect to the ground state that 
we already have. The state Eq. (fT5|) is determined vary- 
ing the energy functional with respect to f^ e and T>\. 
When n 2 — 1 configurations have already been computed 
for the first excited state, one makes the ansatz 



\<p?) = cr*\<P?) 



/3" 2 I02 2 ) 



(16) 



where the state \4>2k) nas a f orm similar to Eq. (fT2"|) but 
written in terms of the new coefficients /J and the GHF- 
transformations T> 2 (i= 1, . . . , n%). Once again, the coef- 
ficients a™ 2 and /3" 2 are obtained by requiring orthonor- 
malization with respect to the ground state we already 
have. The wave function Eq. (|16p is determined vary- 
ing the energy functional with respect to the last added 
GHF-transformation V^ 2 and all the coefficients /J ■ 

Now, we consider the most general situation in which 
the ground state l^" 1 ) = l^™ 1 ) as well as a set of m — 
2 Gram-Schmidt orthonormalized excited states \^ 2 ) 1 
|<^3 3 ), . . . , l^m"!!. 1 ), all of them with the same symmetry 
quantum numbers O, are already at our disposal. Each 
of these m — 1 states is optimized using chains of VAP 
calculations, as discussed in Sec lIIB H and in the present 
section. The key question is then how to approximate 
the m th excited state by n m non-orthogonal symmetry- 
projected configurations. We also need to ensure orthog- 
onality with respect to all the m — 1 states we already 
have. Let us assume that n m — 1 configurations have 
already been computed for the m th excited state with 
symmetry O. Then an approximation in terms of n m 
non-orthogonal symmetry-projected GHF-configurations 
is obtained with the GHF-EXC-FED ansatz 
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FIG. 1: (Color online) The ratio of correlation energies k, 
obtained with the GHF-FED approximation is plotted as a 
function of the inverse of the number of GHF-transformations 
for lattices with N S it es = 20 and 30. Results are shown for 
on-site repulsions of U = 2i, 4i and St. For more details, see 
the main text. 



m— 1 
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\<PrrT) = 2^ W * l"?V)+ T 
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where the state \$%j( ) has a form similar to Eq.(fl2l 
but is written in terms of new coefficients /^f and GHF 
transformations P^ (i= 1, . . . , n m ). The coefficients r" 
and wf read 



Iji-S^ICr)" 



1/2 



W, 



-E(^) tt {Cift ro }r" 

in terms of the projector (i.e., 5 m _i = S^i-i) 



(18) 



where the Hamiltonian •%"' ,,e and norm A/'"™ 6 ker- 
nels account for the fact that m — 1 linearly indepen- 
dent solutions have been removed from the variational 
space. The kernel expressions are slightly more in- 
volved than the ones in Eqs.© and (Tl4l) but still can 
be obtained straightforwardly. The variation with re- 
spect to the coefficients f™ leads to a generalized eigen- 
value equation similar to Eq. (|10l) with the constraint 
/ n ™ e t J\fnm® jn m e _ -^ w hil e the unrestricted minimiza- 
tion with respect to the last added GHF-transformation 
V 7 ^ is carried out via the Thouless theorem. 

The GHF-EXC-FED scheme outlined in this section 
provides, for each set of symmetry quantum numbers 
0, a (truncated) basis of m (orthonormalized) states 
l^™^ 9 ), . . . , l^"^ 3 ), each of them expanded by m, . . . ,n m 
non-orthogonal symmetry-projected GHF-determinants, 
respectively. Finally, the diagonalization of the Hamilto- 
nian Eq.([T]) in such a basis 



m 



C% = Q 



(21) 



j=i 



provides ground and excited states 



in' 



3 = 1 



(22) 



which may account for additional correlations. Never- 
theless, because many of these correlations have already 
been accounted for in the MR expansion of each of the 
m basis states (as discussed above), one may expect the 
role of this final diagonalization to be, in general, less 
important than in the scheme used in Ref. 36. 

Both the GHF-FED and GHF-EXC-FED VAP approx- 
imations could be extended to the use of general Hartrce- 
Fock-Bogoliubov (HFB) transformations. 40 This, how- 
ever, would require an additional particle number sym- 
metry restoration that increases our numerical effort by 
around one order of magnitude and has hence not been 
included in the present study. 



C. Spectral functions and density of states 



5 ro _ 1= J2 IC)(^) tt WI 



.fc=i 



(19) 



T\<PT)- The MR 



with the overlap matrix Aik — Wi Wk 
GHF-EXC-FED wave function Eq.([T7|> is determined by 
varying all the coefficients /^f and the last added GHF- 
transformation 2?"™ . The energy is 



E 



n m B _ J H J 

fn m 0f\fn m fn m & 



(20) 



In this section, we briefly discuss the computation of 
the SFs and DOS within our theoretical framework. Let 
us assume that for an A^-electron system we have already 
obtained, along the lines described in Sec lIIB j"l a GHF- 
FED ground state solution 10™^-" )• For all the lattices 
considered in the present study the ground state has spin 

5 = but not necessarily linear momentum zero [i.e., 

6 = (0, k)]. In all cases, the ground state transforms as 
an irrep of dimension 1. Therefore, for this specific case, 
we can simply write the ground state wave function as 
\ni,N e ,k). The ground state energy will be denoted as 

j^nik 



600 








s^ 


540 
480 




U=4t 


420 




(b) j/ 


' 


o- 360 
£300 






■ 






■ 


240 
180 


Nsites=50 I! 


120 






■ 


60 


'yS^ 


. 


. 



write 



10 15 20 25 









"1 




30 






(a) . 








U=4t 


q.20 

CL 






CO 

10 




N s ites=50 I" 











60 



120 180 240 300 
number of proccessors 



360 



FIG. 2: Speedup of a typical GHF-FED calculation is shown 
in panel (a) as a function of the number of proccessors. The 
corresponding scaling (for a fixed number of proccessors) with 
the number of transformations m is presented in panel (b). 
Results are for Nsites = 50 and U = 4i. 



Usually, the SFs are defined as the imaginary part of 
the time-ordered Green's function and can be calculated 
from the Lehmann representation^ In order to compute 
them, we approximat e 35 ' 36 the ground states of the (iV e ± 
l)-electron systems with the quantum numbers = 
(1/2, fc ± ). For the (N e — l)-electron system we superpose 
the Goldstone (hole) manifolds R(Q)R{g)b h (V\)\V\) in 
the ansatz 



\riT, N e + 1, ft 



+\ - 






PfMbUvl)\Vl) (24) 



ipM 



where the index i runs again as i = l,...,nr and 
p = N e + 1, . . . , 2N S ites- The mixing coefficients f lB 
and g l ® are determined by solving eigenvalue equations 
similar to Eq. tfTT))) . This yields a maximun of 2tlt xN e xd 
hole solutions with energies E nTk and a maximum of 
2tit x (2 A? sites — N e ) x d particle solutions with ener- 
gies E nTk for each irreducible representation of the 
space group. The quantity d is the dimension of the 
corresponding irreducible representations, i.e., d = 1 for 
A;* = 0, 7T and d = 2 for /c ± ^ 0, tt. The hole B(q, to) and 
particle A(q,u)) SFs are then written in their standard 
form 



B(c 7 ,a;)=^KnT I ^e-l,A;"|c ?CT K,iV e ,fc)| 2 

k~ <t 

xS (u - E nik + E nTk ~\ 

fc+(T 

x5 (u - E nTk+ + E nik ) (25) 

and the DOS can be computed as 



tf(u)=J2V 3 (q,u)+A(q,u)) 



(26) 



Due to the finite size of the considered lattices, both 
the hole and particle SFs consist of a finite number of <5 
functions with different weights. Therefore, we introduce 
an artificial Lorentzian width T for each state. 



\n T , N e - 1, ft") = Y, PhM P%Mh(P[)\V[) (23) 

ihM 

where i = 1 , . . . , tit ■ The number ut of GHF- 
transformations used to expand the state Eq. (|23p may 
be different from the one (i.e., n\) in the GHF-FED 
ground state wave function. We write bh(T>\) to ex- 
plicitly indicate that holes are made on different in- 
trinsic determinants \T>\) corresponding to the lowest- 
energy states of the iV e -electron system approximated 
by a single symmetry-projected configuration along the 
lines described in Sec III Al The hole index h runs as 
h = 1, . . . , N e . Note, that the label b = ±1 is not explic- 
itly written in this section, but it is taken into account 
whenever needed. 

For the (7V e + l)-electron system we superpose the 
Goldstone (particle) manifolds R{n)R(g)bl(V[)\V\} and 



III. DISCUSSION OF RESULTS 

In this section, we discuss the results of our calcula- 
tions for some illustrative examples. In most cases, we 
consider on-site repulsions of U = 2t, At and 8t represent- 
ing weak, intermediate-to-strong (i.e., non- interacting 
band width) and strong correlation regimes. First, in 
Sec lIII Al we consider the ground states of half- filled lat- 
tices of various sizes. We compare our ground state and 
correlation energies with the exact ones, as well as with 
those obtained using other theoretical approaches. We 
then discuss the dependence of the predicted correlation 
energies on the number n\ of non-orthogonal symmetry- 
projected configurations used to expand our ground state 
wave functions. The computational performance of our 
scheme is also addressed. The structure of the intrinsic 
GHF-determinants resulting from our GHF-FED YAP 



Z 




FIG. 3: (Color online) The quantity £J [Eg.(|28p] is plotted as a function of lattice site j for some typical symmetry-broken 
GHF-determinants \D\) resulting from the GHF-FED VAP optimization for N attea = 50 and U = At. The UHF spin-density 
wave is plotted in red for comparison. For more details, see the main text. 



procedure is discussed in Sec lIIIBl Our results for SS- 
CFs in real space and MSFs, for lattices with up to 30 
sites, are presented in Sec lIIICI In Sec lIIID( we compare 
the DOS provided by our theoretical framework with the 
exact one, obtained with an in-house full diagonaliza- 
tion code in a lattice with N S i tes = 10. Hole SFs are 
also discussed in the case of N S i tes = 30. Finally, in 
Sec lIIIE| we present results obtained for the excitation 
spectra in lattices with N S i tes — 12, 14 and 20 and also 
discuss the structure of the underlying symmetry-broken 
GHF-determinants resulting from our GHF-EXC-FED 
VAP procedure for excited states. 



A. Ground state and correlation energies 

Let us start by considering lattices of 12 and 20 sites 
with the GHF-FED scheme discussed in Sec lIIBll The 
corresponding O = (0, n) ground states have B\ symme- 
try, i.e., they are symmetric under the reflection x —> — x. 
In Table I, we compare the predicted ground state ener- 
gies with the exact ones. For completeness, we also in- 
clude energies provided by the standard (i.e., one trans- 
formation) restricted (RHF) and unrestricted (UHF) HF 
frameworks. Ours is a VAP approach whose quality can 
be checked by studying how well it reproduces the ex- 
act ground state correlation energies. To this end, we 
consider the ratio 



kghf-fed — 



Erhf — Eghf-fed 
Erhf — Eexact 



x 100 



(27) 



between the GHF-FED and the exact correlation ener- 
gies. For the UHF approximation, kuhf is obtained from 
a similar expression. 

We observe from Table I that the inclusion of ni=10 
non-orthogonal symmetry-projected configurations with 
the GHF-FED approach significantly improves correla- 
tion energies with respect to UHF. In fact, kghf-fed > 
99.79% in all considered correlation regimes even for 



Nsites = 20 which is out of reach with exact diagonal- 
ization. 

In the case of N si tes — 14, whose O = (0, 0) ground 
state has A\ symmetry, i.e., it is symmetric under the 
reflection x — >• —x, our calculations with m = 10 
transformations predict energies of — 11.9539*, —8.0874*, 
and —4.6127* compared to the exact ones of — 11.9543*, 
-8.0883*, and -4.6131* for U = 2t,4t, and 8*, respec- 
tively. This yields kghf-fed values of 99.95, 99.97, and 
99.99%, respectively. Results for this lattice have been 
reported in the literature with the ResHF framework us- 
ing the half-projection method/^ 5 - For half-filled lattices 
with sizes comparable to the ones already mentioned, the 
Gutzwiller method 56 provides k ratios around 85, 77, and 
50%, respectivel y 58 ' 59 Let us also mention that our GHF- 
FED energies for N S i tes = 12 and 14 improve upon pre- 
viously reported VAP values of -6.9093* and -8.0577* 
for U = 4i^ 

Calculations have also been carried out for N si t e s = 16, 
whose O = (0, 7r) ground state has B\ symmetry. We 
have obtained ground state energies of — 16.4754* for 
U = t and —9.2122* for U = 4* while the exact ones 
are — 16.4758* and —9.2144*, respectively. The energies 
reported with the density matrix renormalization group 
(DMRG) method 17 in real space [momentum space ] are 
-16.4755* and -9.2144* [-16.4734* and -8.8349*], re- 
spectively. 

The ground state energies obtained for lattices with 
Nsites — 30 and 50 are compared in Table II with the 
exact ones. In this case, the corresponding O = (0, 0) 
ground states have A\ symmetry. In the same table, 
we also present ground state energies predicted with the 
ResHF method 4 - based on n\ — 30 UHF-transformations 
(i.e., UHF- ResHF). It is very satisfying to observe that 
both the GHF-FED and the UHF-ResHF VAP schemes 
can account for k > 98% in a relatively large lattice with 
Nsites = 30. In fact, the GHF-FED scheme provides 
k > 99.39% with 45,048 variational parameters that rep- 
resents a small fraction of the dimension of the restricted 
(i.e., accounting for all symmetries) Hilbert space in this 
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FIG. 4: (Color online) Ground state spin-spin correlation 
functions in real space for lattices of different sizes. Results 
are shown for \J=2t (a), At (b) and 8t (c). In each panel, 
the inset displays a close-up of the long-range behavior pre- 
dicted by the GHF-FED scheme (red curve) compared with 
the one obtained within the standard UHF approximation 
(green curve) in the case of the N s u es =S0 lattice. For more 
details, see the main text. 



lattice. In this case, our GHF-FED energy also improves 
the variational value — 16.6060£ obtained in Rcf. 35 for 
U = At using a single symmetry-projected configuration. 
Note, that the ResHF metho d 41 ! 42 is not intrinsically lim- 
ited to the use of UHF-transformations and, therefore, 
the UHF-ResHF ground state energies shown in Table 
II can still be improved by, for example, adopting GHF- 



transformations as basic building blocks. 60 

Let us now comment on our results for N S it es — 
32 whose = (0, n) ground state has B\ symme- 
try. We have used n\ = 25 GHF-transformations. 
For on-site repulsions of U = t and 2£, we have ob- 
tained energies of — 33.2137£ and — 26.9814£ while the 
UHF-ResHF ones 4 ^ are -33.2128* and -26.9556t, re- 
spectively. These energies, should be compared with 
the exact ones of — 33.2152i and — 27.0183t as well as 
with the values — 33.2008t and — 26.8016£ obtained with 
the momentum-space DMRG method^- Obviously, the 
previous momentum-space DMRG energies could also 
be improved since a given truncation error depends on 
the number of eigenstates kept in the renormalization 
procedure^ 

From the previous results, we conclude that the GHF- 
FED approximation can be considered a reasonable start- 
ing point for building correlated ground state wave func- 
tions which, at the same time, respect the original sym- 
metries of the ID Hubbard Hamiltonian. This is fur- 
ther corroborated from the results, shown in Table II, 
for N s ites = 50. In particular, even when our descrip- 
tion of the ground state in this lattice is poorer than 
in the N S ites — 30 case since we have kept the same 
number of GHF-transformations, it is remarkable that we 
obtain (with 125,048 variational parameters) the values 
kghf-fed = 96.46,95.19, and 98.65%, respectively. For 
the same on-site repulsions the variational Monte Carlo 
method- 9 predicts k values of around 87, 92, and 96%. 
The corresponding UHF-ResHF values 4 ^ are also listed 
in Table II. 

In Fig[TJ we have plotted the ratio kghf-fed, as a 
function of the inverse 1/ni of the number of transfor- 
mations rii included in the GHF-FED ansatz, for lat- 
tices with N S ites = 20 and 30. They increase smoothly 
with the number of non-orthogonal symmetry-projected 
configurations used to expand the wave function. From 
FigHJ it is apparent (see also Tables I and II) that with 
increasing lattice size, we need a larger number ni of 
symmetry-projected configurations to keep and/or im- 
prove the quality of the GHF-FED wave functions. For 
example, comparing the N S it es = 20 and 30 lattices, 
we see that in the former n\ = 10 transformations are 
enough to obtain kghf-fed > 99.79% while in the lat- 
ter 98.69% < kghf-fed < 99.49%. On the other hand, 
in the N S i tes = 50 case, n\ = 10 transformations leads to 
93.37% < kghf-fed < 97.68%, whereas with m = 25 
we reach the kghf-fed values shown in Table II. 

Obviously, as in many other approaches to many- 
fermion systems, we are always limited to a finite num- 
ber of configurations in practical calculations. Never- 
theless, the GHF-FED scheme provides compact ground 
state wave functions whose quality can be systematically 
improved by adding new (variationally determined) non- 
orthogonal symmetry-projected configurations. In fact, 
both ours and the ResH F 42 i 44 i 45 wave functions are noth- 
ing else than a discretized form of the exact coherent- 
state representation of a fermion stated 2 , and, therefore, 



10 




2.4 2.6 2 



InNc, 



FIG. 5: (Color online) The GHF-FED magnetic structure 
factor is plotted as a function of In Nsites for different lattice 
sizes. GHF-FED results are shown for on-site repulsions of 
U — It (red diamonds), 4i (black diamonds), and 8£ (blue 
diamonds). A straight line has been fitted to guide the eye. 
The corresponding magnetic structure factors predicted by 
the UHF approximation for U — 2t (continuous red curve), 
it (continuous black curve), and 8t (continuous blue curve) 
are also included for comparison purposes. 



become exact in the limit ri\ — > oo. Our aim in the 
present work is not to lower the ground state energy as 
much as possible but to test to which extent our scheme 
can account for relevant correlations in the considered 
lattices. Therefore, for the largest lattices here studied 
(i.e., N S ites — 30 and 50), we have restricted ourselves 
in practice to a maximun number n\ = 25 of GHF- 
transformations. 

A few words concerning the computational perfor- 
mance of our method arc in order here. In panel (a) 
of Figl^l we have plotted the speedup of a typical calcu- 
lation as a function of the number of proccessors. Results 
are shown for N s ues = 50 and U = 4t but similar behav- 
ior was also obtained for U — 2t and 8t. As demonstrated 
in the plot, the GHF-FED speedup grows linearly with 
the number of processors used in the calculations. On the 
other hand, panel (b) of the same figure shows that (for a 
fixed number of proccessors) an efficient implementation 
of our variational scheme scales linearly with the number 
rii of GHF-transformations used while the ResHF scaling 
is quadratic. 



B. Structure of the intrinsic determinants and 
basic units of quantum fluctuations 

An interesting issue is whether there is any relevant 
information in the symmetry-broken (i.e., intrinsic) de- 
terminants \V\) resulting from the GHF-FED VAP opti- 



mization. We are interested in comparing the structure 
of these determinants with the spin-density wave solution 
obtained with the standard UHF approximation. Here, 
one should keep in mind that a variationally optimized 
GHF-determinant has the same energy as the optimal 
UHF one. 63 We have studied the quantity 



CiU) = (-r^iiso'M) • (v{\s(i)\v\) (28) 

where j = l,...,N S it es is the lattice index while i = 
1, . . . , rii enumerates the GHF-determinants in the GHF- 
FED ground state (subscript 1) solution. Among the 
rii = 25 transformations T>\ used for the N S ites~50 lat- 
tice, we have selected some typical examples to plot the 
quantity £\(j). Results are displayed in panels (a) to (d) 
of Fig |3J Other determinants \T>\), not shown in the fig- 
ure, exhibit the same qualitative features. Similar results 
are also found for other on-site repulsions, as well as for 
other lattices. 

For the standard UHF spin-density wave solution, the 
quantity £\(J) has nearly constant positive values plot- 
ted with red lines in FigJ3] A very different behavior ap- 
pears in the intrinsic GHF-determinants \T>\) associated 
with the GHF-FED solution. First, we observe a broad 
spin feature distributed all over the lattice, which is a 
consequence of the richer spin textures provided by the 
use of GHF-transformations and full 3D spin projection. 
In addition, pairs of points (black squares) appear where 
£i(j) changes its sign (i.e., the spin-density wave reverses 
it phase). These defects of the spin-density wave phase 
represent soliton-antisoliton (S — S) pairs in the case of 
half-filled lattices i 42 ' 43 ' 45 ' 52 In particular, our analysis of 
the charge densities p\(j) = 1 — S CT (^ ) il"io-|^ , i) reveal 
that they correspond to neutral S° — S° pairs. Let us 
stress that the presence of at least one S° — S° pair is a 
genuine VAP effect appearing even if we approximate a 
given ground state within a SR framework^ 53 as discussed 
in Seeing 

Furthermore, Figj3] illustrates how the S° — S a pairs 
appear at different lattice locations j with varying dis- 



tance R s o__~qo among the members of the pairs. The 

latter represents the breathing motion of the S° — S a 
pairs. The S° — S° pairs are present in all the intrinsic 
determinants \T>\) associated with the GHF-FED expan- 
sion, which as already mentioned above, superposes the 
Goldstone manifolds Rs(Q)R(g)\'D\) containing defects 
in the spin-density wave. We are then left with an in- 
tuitive physical picture in which the soliton pairs can be 
regarded as basic units of quantum fluctuations in our 
GHF-FED states. On the other hand, the interference 
between S° — S° pairs belonging to different symmetry- 
broken determinants \D\) is accounted for in our calcula- 
tions through a resonon-like equation similar to Eq. (1101) . 
This interpretation has been suggested in previous stud- 
ies with the ResHF method^ 2 - - — 
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FIG. 6: (Color online) The DOS (black) for N a ites — 10 at U — 2t, At, and 8t is plotted in panels (a), (b), and (c), respectively, 
as a function of the shifted excitation energy uj — U/2 (in t units). Results have been obtained by approximating the N e and 
(iV e ± l)-systems with n\ = 10 and ut — 10 GHF-determinants. Our results are hardly distinguishable from the DOS obtained 
with exact diagonalization (red). A Lorentzian folding of width V = 0.05i has been used. For more details, see the main text. 



C. Spin-spin correlation functions and magnetic 
structure factors 

Let us now consider the ground state spin-spin corre- 
lation functions (SSCFs) in real space. For a given set of 
symmetry quantum numbers G, they can be computed 

as 



■^ro') 



6^|S(j).S(l)|^ e ) 

^Tk\4>Tk & ) 



(29) 



Note that if a wave function has good spin, as it is 
the case with the GHF-FED one, the SSCFs have to be 
the same for all the members of a (2S+l)-multiplet and, 
therefore, they cannot depend on the S quantum number. 
However, a dependence with respect to the particular row 
m of the space group irreducible representation that we 
are using in the projection still remains and is explicitly 
included in J"" ie (j). 

The SSCFs corresponding to the ground states for 
N S ites = 14, 18, 22, 26, and 30, approximated by n\ = 
10, 10, 15, 25, and 25 GHF-transformations, respectively, 
are depicted in panels (a), (b) and (c) of FigH) We 
have not considered the N sites — 50 lattice because 
our description of its ground state is poorer than for 
N sltes = 30. The GHF-FED SSCFs display a rapid de- 
crease for j < 3. A similar feature has been studied in 
previous works £2 r 45 M Regardless of the on-site interac- 
tion, the short range part of the SSCFs runs parallel for 
the lattices considered in Fig|4] pointing to converged be- 
havior as a function of lattice size. Morevover, the mid 
and long range amplitude of the SSCF for a given lattice 
increases with increasing U values. 

In each panel of Fig HI the inset displays a close-up 
of the long range behavior of the SSCF predicted by the 
GHF-FED approximation (red curve) compared with the 
one obtained within the standard UHF approach (green 



curve) for N sit es = 30. As can be observed, the am- 
plitude of the UHF SSCF remains constant for j > 5 
while the GHF-FED one exhibits a damped long range 
trend. Previous studies have suggested that there are two 
important ingredients necessary to account for a quali- 
tatively correct long range behavior of the SSCFs: the 
self-consistent optimization of the intrinsic determinants 
(i.e., orbital relaxatio n 43 i 45 ) and having pure spin states 
(i.e., no spin contamination 4 ^). Our wave functions meet 
both conditions. 

The magnetic structure factors (MSFs) can be com- 
puted as 



™i6 



N s 



ij 



R 



i+j _ 



.-fis^-sm^f) 



-me 1 

J 1K 



ini0 



(30) 



and the ones corresponding to the ground states for 
N sltes = 14, 18, 22, 26, and 30 are displayed in Figf5] as 
functions of In N S i tes . The UHF MSFs are also included 
in the plot for comparison purposes. At variance with the 
UHF MSFs which diverge exponentially, our GHF-FED 
results display an almost linear behavior. A previous 
work 65 has shown that the SSCFs in real space behave 
for a half- filled system as ~ (ln° j)/ j. This implies that 
as a funcion of the lattice size, the MSFs should behave 
as ln 1+rT N S ites- In FigfSJ we have simply fitted a straight 
line to guide the eye. We have not attempted to deter- 
mine logarithmic corrections as this would require larger 
lattices than those studied in the present paper. 



D. Spectral functions and density of states 

In panels (a), (b), and (c) of FigEl we have plotted 
(black) the DOS J\f(u) for N S ites = 10. In the same fig- 
ure, we have also plotted (red) the exact DOS obtained 
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FIG. 7: (Color online) The hole SFs for N a ites — 30 with U — 2£, At, and 8i are plotted in panels (a), (b), and (c) as functions of 
the shifted excitation energy w — U/2 (in i units). Results have been obtained by approximating the N B and (N e ± l)-systems 
with m — 25 and nr — 25 GHF-determinants. The hole SFs for momenta identical to the Fermi momentum are displayed in 
brown color. The shapes of some selected hole SFs (i.e., k — 0, 27r/5 and 7n/15), obtained by approximating the ground states 
of the (iV e ± l)-systems with nr = 5 (blue), 15 (red) and 25 (black) but the ground state of the JVe-system always with ni=25 
GHF-transformations, are compared in panels (d), (e), and (f). A Lorentzian folding of width F = 0.05£ has been used. For 
more details, see the main text. 



with an in-house full diagonalization code. There is excel- 
lent agreement between ours and the exact DOS concern- 
ing the position and relative heights of all the prominent 
peaks. Both ours and the exact DOS exhibit the particle- 
hole symmetry well known for half-filled systems^ and a 
splitting into the lower and upper Hubbard bands. The 
Hubbard gap between these bands increases with larger 
U. Both dynamical cluster approximation^— and cel- 
lular dynamical mean field theorySS studies suggest that 
this gap is preserved for any finite value of the on-site 
interaction at sufficiently low temperatures even in the 
thermodynamic limit, with U = Ot being the only singu- 
lar point. 

Tendencies to spin-charge separation as well as other 
relevant shadow features inside the Brillouin zone, similar 
to the ones expected in the infinite-t/ limit of the ID 
Hubbard model^ have been found in previous cellular 
dynamical mean field theory^ and cluster perturbation 
theory 6 ^ studies of the spectral weigths in the case of 
finite on-site repulsions. In panels (a), (b), and (c) of 
Fig 13 we have plotted the hole SFs for the N sit es = 30 
lattice. 

The first feature observed from Fig(3 is the Hubbard 
gap opening at the Fermi momentum kp — 7tt/15. The 
spectral weight concentrates on the prominent peaks be- 
longing to the spinon band. Our calculations for smaller 
lattices with N S i tes = 14 and 20 indicate that this spinon 
band is quite stable in terms of lattice size although 
the relative height of its peaks decreases for increas- 



ing U values. The holon singularities are clearly visible 
in some of the SFs shown in Figj3 for linear momenta 
— 7r/2 < k < tt/2. On the other hand, the holon bands 
can also be followed for linear momenta k > tt/2 and 
k < —tt/2. They are the mirror images of the ones 
with opposite w — U/2 values^ and become apparent 
for U = At and 8i. However, besides the spinon band, 
the most relevant feature in our SFs is the very extended 
distribution of the spectral weight for linear momenta 
—w/2 < k < tt/2. The comparison with our SFs for 
N S ites = 14 and 20, obtained with n\ = 10 and nr = 10 
GHF-transformations, reveals that the increase of lattice 
size produces more pronounced shadow features due to 
the fragmentation of the spectral strength over a wider 
interval of to — U/2 values. The previous finite size re- 
sults show that our SFs exhibit tendencies beyond a sim- 
ple quasiparticle distribution and agree qualitatively well 
with the ones obtained using other approximations! 66 ' 68 
The shapes of some selected hole SFs are compared in 
panels (d), (e), and (f) of Fig 13 As can be seen from 
panel (d), nx = 15 transformations are enough to ac- 
count for all relevant details of the SFs shown in panel 
(a) for U — 2t. On the other hand, panels (e) and (f) 
show that a larger number of transformations is required 
for U = At and 8t. In particular, increasing the number 
of transformations for the (N e ± l)-systems from ny = 5 
to ut = 15 and/or 25 leads to a shift of the main peaks 
and redistributes the spectral strength of some of the 
peaks found in the SFs for ut = 5 as a result of the 
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FIG. 8: The energies of some selected states obtained within 
the GHF-EXC-FED approximation for N sites = 12 and 14 
are compared with the ones provided by single-reference VAP 
(SR VAP) calculations (only 3D spin and linear momentum 
projections) as well as with the exact ones from Lanczos 
diagonalization. 35 Results are shown for U — At. 



small number of configurations used in the calculations. 
This explains the differences between ours and the SFs 
reported, for the same lattice at U — At, in the previous 
VAP study2£ of the ID Hubbard model using only n\ = 1 
and tit = 5 GHF-transformations. 



E. Excitation spectra 

In this section, we consider the low-lying excitation 
spectra obtained for N 'sites = 12, 14, and 20 with the 
GHF-EXC-FED scheme discussed in Sec lIIB21 For each 
irreducible representation of the space group, we have 
computed the lowest-energy and first excited states with 
spins S — 0,1, and 2. Each of these states has been 
approximated by 10 non-orthogonal symmetry-projected 
GHF-detcrminants. A final 2x2 diagonalization of the 
ID Hubbard Hamiltonian has also been carried out. For 
these particular lattices, we have found that for each sym- 
metry 0, the (Gram-Schmidt orthonormalized) ground 
|y™^-~~ '"} and first excited \^> 2 k~ ' ) states are very 
weakly coupled through the Hamiltonian. Due to this, 
the energies corresponding to the states \&>f K ) and |Slf A -) 
resulting from the 2x2 diagonalization are almost identi- 
cal to those corresponding to the basis states \<Pi^~ ' ) 
and \f2K~ ' )■ However, this cannot be anticipated a 
priori and the final diagonalization of the Hamiltonian 
should always be carried out. 

In FigJHl we compare the energies of some selected 



states for N s u e s = 12 and 14 with the ones obtained 
in the previous variational study^ of the ID Hubbard 
model where 3D spin and linear momentum projections 
were carried out. The exact results in Fig|S] correspond 
to Lanczos diagonalizations. As can be observed, our MR 
calculations, where in addition to 3D spin projection the 
full space group of the ID Hubbard model is taken into 
account, improve the energies reported in Ref. 35 for 
both ground and excited states. 

In Fig[31 we show the low-lying spectrum obtained 
via Eq. ([2T|) . for N s ues — 14 [panels (a), (b), and (c)] 
and 20 [panels (d), (e), and (f)] taken as representatives 
examples of systems whose 9 = (0,0) and 6 = (0,7r) 
ground states have A\ and B\ symmetries, respectively. 
We observe that both the lowest-lying singlet and triplet 
states obtained in our calculations nicely follow the sine- 
like dispersion trend in the exact curve for N S ites —> 
oo. The anomaly observed in the GHF-EXC-FED k- 
dispersion for the lowest-energy singlets and triplets has 
also been found in previous studies within the ResHF 
framework^ as well as in finite versions of the exact 
Licb-Wu solutions^ For any finite U value, the exact 
Nsites -^ oo curves exhibit gapped excitations, excep- 
tion made for the k = and k — n states which are 
degenerate. In our calculations such a degeneracy is bro- 
ken due to finite size effects. However, we observe that 
for a given finite U value, the energy difference between 
the lowest singlet and triplet states decreases with in- 
creasing lattice size. For example, for U = 2t we have 
obtained AE s -t — 0.5287t in the case N S ites — 14 while 
AE s -t — 0.1275t for N s ues — 20. For increasing on- 
site respulsions, irrespective of the lattice size, an overall 
compression of the spectra takes place. This is consistent 
with the fact that in the limit U — >■ oo all the configura- 
tions shown in Figj9] should become degenerate. 



F. Structure of the intrinsic determinants and 

basic units of quantum fluctuations in the 

GHF-EXC-FED wave functions 



In Sec lIIIBl we have discussed the structure of the 
intrinsic determinants associated with the GHF-FED 
states. Here, we pay attention to the symmetry-broken 
ones used to expand the GHF-EXC-FED wave functions. 
To illustrate our results, we consider states belonging to 
the spectrum shown in panel (e) of FigJSJ In particular, 
we have plotted in panels (a) to (d) of FigfTU]the quanti- 
ties £x(i) (black) and QU) ( rec computed [see, Eq.(|2"5])] 
with some of the n\ — 10 and ri2 = 10 symmetry-broken 
determinants \T>\) and \T> 2 ) used to expand the lowest- 
energy Ypxk ' ~ ) and first excited Yp^k ' ~ ) states with 
© = (1,0) and Ai symmetry. Other determinants, not 
shown in the figure, exhibit the same qualitative features. 
Similar results are also obtained for U = 2t and 8t as well 
as for other lattices. 

As can be observed, both Q(j) and £|(j) display de- 
fects similar to the ones already discussed for the S = 
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FIG. 9: (Color online) Energy spectrum obtained via Eq. (|2ip for N 3 ites — 14 [panels (a), (b), and (c)] and 20 [panels (d), (e), 
and (f)]. For each irreducible representation of the space group, the lowest-energy and first excited states with the spins S — 
(red bars),l (blue bars), and 2 (green bars) have been plotted. Results are shown for U = 2t,4t, and St. The exact dispersion 
curves^ for N S it es — > oo (thin black lines) are also included. In order to guide the eye, the lowest-lying states with spin 5 = 
(S = 1) have been connected by long (short) dashed lines. For more details, see the main text. 



ground states provided by the GHF-FED approximation 
(see, Fig|3]). From this we conclude that not only the 
ground but also the excited state wave functions provided 
by our MR VAP scheme superpose Goldstone manifolds 
built in terms of intrinsic GHF-determinants containing 
defects (i.e., solitons) that can be regarded as basic units 
of quantum fluctuations. In general, the intrinsic de- 
terminants associated with different symmetry-projected 
states may develop local variations of the charge density 
as seen (dashed blue curve) from Fig llOl where we have 
also plotted the quantity p\(j) = 1 — X^^l^jcrl^-'l)- 



IV. CONCLUSIONS 

The accurate description of the most relevant correla- 
tions in the ground and low- lying excited states of a given 
many-fermion system, with as few configurations as pos- 
sible, is a central problem in quantum chemistry, solid 
state, and nuclear structure physics. In the present study, 
we have explored a VAP MR avenue for the ID Hubbard 
model. The main acomplishments of the present work 
are listed below. 

(i) We have presented a powerful methodology of a 
VAP MR configuration mixing scheme, originally devised 
for the nuclear many-body problem, but not yet consid- 



ered to study ground and excited states, with well defined 
symmetry quantum numbers, of the ID Hubbard model 
with nearest-neighbor hopping and periodic boundary 
conditions. Both ground and excited states are expanded 
in terms of non-orthogonal and Ritz-variationally opti- 
mized symmetry-projected configurations. The simple 
structure of our projected states allows an efficient par- 
allelization of our variational scheme, which scales lin- 
early with the number of processors as well as with the 
number of transformations used in the calculations. The 
method also provides a (truncated) basis consisting of 
a few Gram-Schmidt orthonormalized states. This ba- 
sis may be used to diagonalized the Hamiltonian to ac- 
count, in a similar fashion, for additional correlations in 
the ground and excited states with well defined symmetry 
quantum numbers. 

(ii) We have shown that our MR approximation gives 
accurate ground state energies and correlation energies 
as compared with the exact Lieb-Wu solutions for rel- 
atively large half-filled lattices up to 30 and 50 sites. 
The comparison with other theoretical approaches also 
reveals that our scheme can be considered as a reason- 
able starting point for obtaining correlated ground state 
wave functions in the case of the ID Hubbard model. 
We have computed the full low-lying spectrum for the 
Nsites = 14 and 20 lattices. The momentum dispersion 
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FIG. 10: (Color online) Structure of some typical symmetry-broken GHF-determinants \D\) (black) and \T>\) (red) used to 
expand the lowest-energy and first excited states with spin S — 1, linear momentum k = and Ai symmetry. The charge 
density (dashed blue) corresponding to the determinants \T>2) is also included in each panel. Results are shown for U = 4£. For 
more details, see the main text. 



of the lowest-lying singlet and triplet states follows the 
exact shape predicted by the Lieb-Wu solution in the 
thermodynamic limit. With increasing U we also observe 
a general compresion of the spectrum. 

iii) From the analysis of the structure of the intrinsic 
determinants associated with our MR ground and excited 
state wave functions, we observe that they all contain 
defects (i.e., solitons) that can be regarded as basic units 
of quantum fluctuations for the considered lattices. 

(iv) Our results for the ground state SSCFs in real 
space show long range decay that is not observed in the 
UHF case. The MSFs computed from such correlation 
functions show a behavior approximately linear in In 
N s ites consistent with previous results available in the 
literature. 

(v) Our approximation also allows to compute SFs and 
the DOS. To this end, we considered ansatze, whose flex- 
ibility is determined by the numbers n\ and nr of HF- 
transformations used to expand the wave functions of 
systems with N e and (N e ± l)-electrons. For a small 
lattice with N S i tes — 10 we have compared the DOS pre- 
dicted within our approach with the one obtained using 
a full diagonalization and found an excellent agreement 
between both. For a larger lattice with N S i tes — 30 our 
scheme provides hole SFs that agree qualitatively well 
with the ones obtained with other approximations and 
exhibit tendencies beyond a simple quasiparticle distri- 
bution. 

We believe that the finite size calculations discussed 
in the present study already show that VAP approx- 
imations, based on MR expansions in terms of non- 
orthogonal symmetry-projected HF-determinants, rep- 
resent useful tools that complement other existing ap- 
proaches to study the physics of low-dimensional corre- 
lated electronic systems. Within this context, the scheme 
presented in this work leaves ample space for further 
improvements and research. First, the number of non- 



orthogonal symmetry-projected configurations used in 
the corresponding MR expansions can be increased to 
improve the quality of our wave functions. Second, we 
could still incorporate particle number symmetry break- 
ing (i.e., general HFB-transformations) and restoration 
(i.e., particle number projection) to access even more cor- 
relations. Third, our scheme can be easily extended to 
the 2D case as well as to doped systems with arbitrary 
on-site interaction strengths. Our approximation is also 
general enough so as to be implemented for the molec- 
ular Hamiltonian 71 as well as for lattices like the hon- 
eycomb, the Kagome or the Shastry-SutherlandtZ^ ones. 
The same VAP MR scheme can also be applied to study 
frustrated Hubbard models in the ID and 2D cases. Fi- 
nally, the MR scheme discussed in the present work could 
also be used as a powerful solver in the framework of 
fragment-bath embedding approximations.- 7 In particu- 
lar, it could replace exact diagonalizations for fragment 
sizes where it is not feasible while still providing highly 
correlated (fragment) wave functions. Obviously, a care- 
ful analysis of the corresponding symmetries should be 
carried out in each case. Work along these avenues is in 
progress and will be reported elsewhere. 
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